| Teamleden | Kaggle Username | GitHub Username |
|---|---|---|
| Busse Heemskerk | bussejheemskerk | BJHeemskerk |
| Declan van den Hoek | declanvdh | DeclanvandenHoek |
| Isa Dijkstra | isadijkstra | IsaD01 |
In dit notebook worden er, aan de hand van een train dataset, verschillende Machine Learning modellen opgesteld om de meest accurate manier te vinden om de verhuuraantallen te voorspellen. Er is niet bekend gemaakt wat het product is dat wordt verhuurd. Voor het project hebben we gewerkt in GitHub, om makkelijk de bestanden te delen. Van elk model zijn de voorspellingen ook geupload naar Kaggle.
# Standaard libaries voor visualiseren en dataverwerking
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
# Libaries gebruikt bij de seizoens analyse
from statsmodels.tsa.seasonal import seasonal_decompose
from scipy.signal import find_peaks
from statsmodels.tsa.deterministic import CalendarFourier
from statsmodels.tsa.deterministic import DeterministicProcess
# Libaries gebruikt bij trend analyse
# Libaries gebruikt bij autocorrelatie
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# Libaries gebruikt bij functies rondom modelleren
import datetime
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import train_test_split
from sklearn.model_selection import TimeSeriesSplit
from statsmodels.tsa.stattools import adfuller
# Libaries van de modellen
from sklearn.linear_model import LinearRegression
from sklearn.neighbors import KNeighborsRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.ensemble import StackingRegressor
from xgboost import XGBRegressor
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.statespace.sarimax import SARIMAX
from prophet import Prophet
Aan de hand van de pandas libary is het nu mogelijk om de csv bestanden in te laden.
# Inladen van de test dataset
data_test = pd.read_csv('test.csv', sep = ',')
# Inladen van de train dataset
data_train = pd.read_csv('train.csv', sep = ',')
Nu de data is ingeladen zullen de eerste 5 regels worden getoond van de train dataset.
# Tonen van de eerste 5 regels
data_train.head()
| date_hour | holiday | weathersit | temp | atemp | hum | windspeed | cnt | |
|---|---|---|---|---|---|---|---|---|
| 0 | 2011-01-01 00:00:00 | 0 | 1 | 0.24 | 0.2879 | 0.81 | 0.0 | 16 |
| 1 | 2011-01-01 01:00:00 | 0 | 1 | 0.22 | 0.2727 | 0.80 | 0.0 | 40 |
| 2 | 2011-01-01 02:00:00 | 0 | 1 | 0.22 | 0.2727 | 0.80 | 0.0 | 32 |
| 3 | 2011-01-01 03:00:00 | 0 | 1 | 0.24 | 0.2879 | 0.75 | 0.0 | 13 |
| 4 | 2011-01-01 04:00:00 | 0 | 1 | 0.24 | 0.2879 | 0.75 | 0.0 | 1 |
In dit hoofdstuk zal er worden gewerkt aan de Explorative Data Analysis. Deze stap zal ons helpen om beter inzichten krijgen met betrekking tot de data.
De minimale vereisten voor Machine Learning met Scikit Learn zijn:
Omdat de data is ingelezen via Pandas staat het al in een Pandas DataFrame, hierdoor kan gebruik gemaakt worden van de .info() methode van DataFrames om te kijken naar missende waarden en de datatypen. Vervolgens kan de .describe() method worden gebruikt om een inzicht te krijgen van hoe de numerieke waarden zijn verdeeld.
# Tonen van de kern informatie over kolommen
data_train.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 16637 entries, 0 to 16636 Data columns (total 8 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 date_hour 16637 non-null object 1 holiday 16637 non-null int64 2 weathersit 16637 non-null int64 3 temp 16637 non-null float64 4 atemp 16637 non-null float64 5 hum 16637 non-null float64 6 windspeed 16637 non-null float64 7 cnt 16637 non-null int64 dtypes: float64(4), int64(3), object(1) memory usage: 1.0+ MB
Uit deze diagram is gebleken dat de data_hour kolom niet een data_time object is.
data_train['date_hour'] = pd.to_datetime(data_train['date_hour'])
data_train.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 16637 entries, 0 to 16636 Data columns (total 8 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 date_hour 16637 non-null datetime64[ns] 1 holiday 16637 non-null int64 2 weathersit 16637 non-null int64 3 temp 16637 non-null float64 4 atemp 16637 non-null float64 5 hum 16637 non-null float64 6 windspeed 16637 non-null float64 7 cnt 16637 non-null int64 dtypes: datetime64[ns](1), float64(4), int64(3) memory usage: 1.0 MB
Uit de .info() diagram is ook af te leiden dat elke cel hetzelfde aantal waardes bevat. Dit zou aan kunnen geven dat er geen waardes missen. De datatypes zijn ook allemaal integers of floats (behalve de datum), hieruit is af te leiden dat deze kolommen bestaan uit numerieke waardes.
# Tonen van informatie over numerieke kolommen
data_train.describe()
| date_hour | holiday | weathersit | temp | atemp | hum | windspeed | cnt | |
|---|---|---|---|---|---|---|---|---|
| count | 16637 | 16637.000000 | 16637.000000 | 16637.000000 | 16637.000000 | 16637.000000 | 16637.000000 | 16637.000000 |
| mean | 2011-12-18 02:20:08.294764544 | 0.028671 | 1.415580 | 0.504745 | 0.482608 | 0.624756 | 0.190310 | 190.477009 |
| min | 2011-01-01 00:00:00 | 0.000000 | 1.000000 | 0.020000 | 0.000000 | 0.000000 | 0.000000 | 1.000000 |
| 25% | 2011-06-27 05:00:00 | 0.000000 | 1.000000 | 0.340000 | 0.333300 | 0.470000 | 0.104500 | 41.000000 |
| 50% | 2011-12-18 06:00:00 | 0.000000 | 1.000000 | 0.520000 | 0.500000 | 0.620000 | 0.194000 | 143.000000 |
| 75% | 2012-06-09 02:00:00 | 0.000000 | 2.000000 | 0.660000 | 0.621200 | 0.780000 | 0.253700 | 282.000000 |
| max | 2012-11-30 23:00:00 | 1.000000 | 4.000000 | 1.000000 | 1.000000 | 1.000000 | 0.850700 | 977.000000 |
| std | NaN | 0.166885 | 0.637298 | 0.192369 | 0.171557 | 0.193227 | 0.121915 | 182.026755 |
Volgens de opdrachtbeschrijving zijn sommige kolommen al genormaliseerd, en zouden er geen maximale waardes boven de 1 moeten zijn. Als we naar de dataset kijken klopt dit, voor de kolommen "temp". "atemp" en "hum". Opvallend is dat de kolom "windspeed" ook genormaliseerd zou moeten zijn, maar de maximale waarde onder de 1.0 ligt. Dit kan aan geven dat de data niet goed is genormaliseerd.
min_value = data_train['windspeed'].min()
max_value = data_train['windspeed'].max()
data_train['windspeed'] = (data_train['windspeed'] - min_value) / (max_value - min_value)
data_train.describe()
| date_hour | holiday | weathersit | temp | atemp | hum | windspeed | cnt | |
|---|---|---|---|---|---|---|---|---|
| count | 16637 | 16637.000000 | 16637.000000 | 16637.000000 | 16637.000000 | 16637.000000 | 16637.000000 | 16637.000000 |
| mean | 2011-12-18 02:20:08.294764544 | 0.028671 | 1.415580 | 0.504745 | 0.482608 | 0.624756 | 0.223709 | 190.477009 |
| min | 2011-01-01 00:00:00 | 0.000000 | 1.000000 | 0.020000 | 0.000000 | 0.000000 | 0.000000 | 1.000000 |
| 25% | 2011-06-27 05:00:00 | 0.000000 | 1.000000 | 0.340000 | 0.333300 | 0.470000 | 0.122840 | 41.000000 |
| 50% | 2011-12-18 06:00:00 | 0.000000 | 1.000000 | 0.520000 | 0.500000 | 0.620000 | 0.228047 | 143.000000 |
| 75% | 2012-06-09 02:00:00 | 0.000000 | 2.000000 | 0.660000 | 0.621200 | 0.780000 | 0.298225 | 282.000000 |
| max | 2012-11-30 23:00:00 | 1.000000 | 4.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 977.000000 |
| std | NaN | 0.166885 | 0.637298 | 0.192369 | 0.171557 | 0.193227 | 0.143312 | 182.026755 |
def boxplot_outliers(dataframe, column_name):
plt.figure(figsize=(8, 6))
plt.boxplot(dataframe[column_name], vert=False)
plt.title(column_name)
plt.show()
for col in data_train.columns:
if col != 'date_hour':
boxplot_outliers(data_train, col)
In de boxplot voor windspeed is te zien dat er een aantal outliers zijn. We gaan kijken naar deze rijen.
outliers_windspeed = data_train[data_train['windspeed'] > 0.58]
outliers_windspeed
| date_hour | holiday | weathersit | temp | atemp | hum | windspeed | cnt | |
|---|---|---|---|---|---|---|---|---|
| 178 | 2011-01-08 17:00:00 | 0 | 1 | 0.16 | 0.1212 | 0.37 | 0.649112 | 69 |
| 194 | 2011-01-09 09:00:00 | 0 | 1 | 0.12 | 0.0758 | 0.46 | 0.614083 | 19 |
| 196 | 2011-01-09 11:00:00 | 0 | 1 | 0.16 | 0.1212 | 0.40 | 0.614083 | 49 |
| 265 | 2011-01-12 12:00:00 | 0 | 1 | 0.20 | 0.1515 | 0.47 | 0.684260 | 55 |
| 271 | 2011-01-12 18:00:00 | 0 | 1 | 0.20 | 0.1515 | 0.47 | 0.614083 | 137 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 15876 | 2012-10-28 17:00:00 | 0 | 3 | 0.46 | 0.4545 | 0.77 | 0.614083 | 225 |
| 16207 | 2012-11-13 01:00:00 | 0 | 3 | 0.44 | 0.4394 | 0.88 | 0.754438 | 5 |
| 16468 | 2012-11-23 22:00:00 | 0 | 1 | 0.44 | 0.4394 | 0.33 | 0.614083 | 79 |
| 16472 | 2012-11-24 02:00:00 | 0 | 1 | 0.32 | 0.2727 | 0.39 | 0.719290 | 32 |
| 16482 | 2012-11-24 12:00:00 | 0 | 2 | 0.30 | 0.2576 | 0.36 | 0.719290 | 199 |
217 rows × 8 columns
class Grafieken:
"""
Een class voor eenvoudig visualiseren met Matplotlib
en seaborn.
Parameters:
----------
df : pandas.DataFrame
Een pandas.DataFrame die data bevat om mee te
visualiseren.
"""
def __init__(self, df, target=False):
"""
De constructor voor de Grafieken class
Parameters:
----------
df : pandas.DataFrame
Een pandas.DataFrame die data bevat om mee te
visualiseren.
target : str (default = False)
De kolomnaam van de target kolom, indien toepasbaar
"""
# Aanmaken self.df
self.df = df.drop('date_hour', axis=1)
# Zetten van de target
if target:
self.target = target
# Aanmaken grafiek standaarden
self.figsize = (15,5)
self.pallete = 'deep'
def boxplots(self):
"""
Maakt een aantal box plots gelijk aan het aantal kolommen
in de Dataframe. De box zijn tegen één specifieke kolom
opgezet, om de verdeling t.o.v. die kolom te tonen.
Parameters:
----------------
df : Pandas DataFrame
Een pandas DataFrame met kolommen waarvan je
de verdeling wilt onderzoeken.
kolom : str
De naam van de kolom waarbij de verdeling bekeken
word, ten opzichte van de andere kolommen.
Returns:
---------------
None :
In plaats van iets te returnen laat de
functie meerdere boxplots zien.
"""
# For-loop om elke kolom een eigen boxplot te geven
for col in self.df.columns:
# if-statement om te voorkomen dat de
# kolom-parameter ook een boxplot krijgt
if col != self.target:
# Zetten van de grootte van de plot
plt.subplots(figsize=self.figsize)
# Maken en benoemen van de assen van de boxplots
sns.boxplot(x=col, y=self.target,
data=self.df, color='pink')
plt.xlabel(col)
plt.ylabel(self.target)
plt.title('{} vs. {}'.format(col, self.target))
plt.show()
def barplots(self):
"""
Maakt een aantal bar plots gelijk aan het aantal kolommen
in de Dataframe. De bar is tegen één specifieke kolom
opgezet, om de verdeling t.o.v. die kolom te tonen.
Parameters:
----------------
df : Pandas DataFrame
Een pandas DataFrame met kolommen waarvan
je de verdeling wilt onderzoeken.
kolom : str
De naam van de kolom waarbij de verdeling
bekeken word, ten opzichte van de andere kolommen.
Returns:
---------------
None :
In plaats van iets te returnen laat
de functie meerdere barplots zien.
"""
# For-loop om elke kolom een eigen barplot te geven
for col in self.df.columns:
# if-statement om te voorkomen dat de
# kolom-parameter ook een barplot krijgt
if col != self.target:
# Zetten van de grootte van de plot
plt.subplots(figsize=self.figsize)
# Maken en benoemen van de assen van de barplot
sns.barplot(x=col, y=self.target,
data=self.df, color= 'skyblue')
plt.xlabel(col)
plt.ylabel(self.target)
plt.title('{} vs. {}'.format(col, self.target))
plt.show()
def lineplots(self):
"""
Maakt een aantal line plots gelijk aan het aantal kolommen
in de Dataframe. De line is tegen één specifieke kolom
opgezet, om de verdeling t.o.v. die kolom te tonen.
Parameters:
----------------
df : Pandas DataFrame
Een pandas DataFrame met kolommen waarvan
je de verdeling wilt onderzoeken.
kolom : str
De naam van de kolom waarbij de verdeling
bekeken word, ten opzichte van de andere kolommen.
Returns:
---------------
None :
In plaats van iets te returnen laat
de functie meerdere lineplots zien.
"""
# For-loop om elke kolom een eigen barplot te geven
for col in self.df.columns:
# if-statement om te voorkomen dat de
# kolom-parameter ook een barplot krijgt
if col != self.target:
# Zetten van de grootte van de plot
plt.subplots(figsize=self.figsize)
# Maken en benoemen van de assen van de barplot
sns.lineplot(x=col, y=self.target,
data=self.df, color= 'green')
plt.xlabel(col)
plt.ylabel(self.target)
plt.title('{} vs. {}'.format(col, self.target))
plt.show()
def heatmap(self):
"""
Deze functie maakt een heatmap op basis
van de DataFrame waarmee de class is gemaakt
"""
# Aanmaken Heatmap
fig, ax = plt.subplots(figsize=(5,15))
ax = sns.heatmap(self.df.corr()[[self.target]].sort_values(
by=[self.target], ascending=False
),square=False, annot=True, cmap='RdYlGn')
plt.show()
def timeplot(self, x_kolom, hue=None):
"""
Maakt een bar en een boxplot op basis van een gegeven tijd.
Parameters:
----------
x_kolom : str
De naam van de tijd gebaseerde kolom
hue : str
De naam van de kolom waarop de bars en boxes
worden ingedeeld
Returns:
----------
None :
In plaats van iets te returnen laat
de functie een bar- en een boxplot zien.
"""
# Aanmaken van plotgrootte en barplot
plt.subplots(figsize=self.figsize)
sns.barplot(data=self.df, x=x_kolom, y=self.target,
hue=hue, errorbar=None,
palette=self.pallete)
# Aanmaken van de assen
plt.title(f'Verhuuraantallen per {x_kolom}')
plt.ylabel('Verhuuraantallen')
plt.show()
# Aanmaken van plotgrootte en boxplot
plt.subplots(figsize=self.figsize)
sns.boxplot(data=self.df, x=x_kolom, y=self.target,
hue=hue, palette=self.pallete)
# Aanmaken van de assen
plt.title(f'Verdeling van Verhuuraantallen per {x_kolom}')
plt.ylabel('Verhuuraantallen')
plt.show()
graf = Grafieken(data_train, 'cnt')
graf.heatmap()
Aan de hand van deze heatmap kun je zien dat de hum, temp en atemp de meeste invloed hebben op de hoeveelheid die wordt verhuurd. Door de frequentie van de data (elk uur) is er een te grote schommeling tussen de mogelijke waarden. Om beter te kunnen kijken naar eventuele patronen is het dus nodig om in te gaan zoomen op bepaalde perioden. Om te beginnen met inzoomen kijken we naar de maand april.
graf.boxplots()
graf.barplots()
graf.lineplots()
In de visualisaties kun je zien dat holiday (0) zorgt voor meer verhuuraantallen van het product. Ook zie je dat hoe hoger de weathersite, hoe meer er van het product wordt verhuurd. Ook zie je een stijgende lijn in de temperatuur en de verhuurde aantallen. Windspeed loopt geleidelijk met een paar uitschieters. Hum loopt juist af, dus een lagere hum zorgt voor meer verhuur, maar de hum moet ook weer niet te laag zijn want dan verhuur je bijna niks.
Omdat er met Time-Series data wordt gewerkt, is het nodig om deze verder te gaan onderzoeken. Hiervoor is eerst een nieuw dataframe nodig, waarin alleen de tijd en de verhuuraantallen staan.
# Aanmaken van nieuw DataFrame
df_time = data_train[['date_hour', 'cnt']].copy()
df_time['date_hour'] = pd.to_datetime(df_time['date_hour'])
df_time = df_time.set_index('date_hour')
df_time.head()
| cnt | |
|---|---|
| date_hour | |
| 2011-01-01 00:00:00 | 16 |
| 2011-01-01 01:00:00 | 40 |
| 2011-01-01 02:00:00 | 32 |
| 2011-01-01 03:00:00 | 13 |
| 2011-01-01 04:00:00 | 1 |
Nu het dataframe df_time is aangemaakt kan deze gebruikt worden voor verschillende time-series analyse technieken. Deze technieken worden gebruikt om patronen uit de tijds data te kunnen halen. Om hiermee te beginnen wordt het dataframe in een grafiek gezet. Hierbij zal 2012 gehighlight worden om snel te kunnen zien waar in de grafiek de strat van het nieuwe jaar zich bevindt.
# Aanmaken van standaard plot
ax = df_time.plot(figsize=(15, 5))
# Toevoegen van shading om beter te zien waar 2012 begint
ax.axvspan('2012-01-01', '2013-01-01', color='green', alpha=0.15)
# Aanpassen van de plot labels en titel
ax.set_xlabel("Datum")
ax.set_ylabel("Verhuuraantallen")
ax.set_title("Verhuuraantallen per uur")
# Tonen van de plot
plt.show()
# Aanmaken van standaard plot
ax = df_time.plot(figsize=(15, 5))
# Toevoegen van shading om beter te zien welke data April is
ax.axvspan('2012-04-01', '2012-05-01', color='green', alpha=0.15)
# Aanpassen van de plot labels en titel
ax.set_xlim('2012-03-22', '2012-05-08')
ax.set_xlabel("Datum")
ax.set_ylabel("Verhuuraantallen")
ax.set_title("Verhuuraantallen per uur van 22 maart tot 8 mei")
# Tonen van de plot
plt.show()
Wat hier te zien valt is dat er een patroon vormt in de pieken. Dit patroon herhaalt zich relatief consistent. Hoewel deze rond 22 april zich niet herhaald, is het belangrijk om niet van de uitzondering de regel te gaan maken. Omdat deze herhaling per maand meerdere malen voorkomt is het van belang om te kijken naar de data per week. Op deze manier kan vastgesteld worden of er per week ook een herlaing plaats vind. Om dit gemakkelijk te berekenen en te visualiseren wordt er gebruik gemaakt van de seasonal_decompose functie van de statsmodels libary.
# Gebruiken van seasonal_decompose om seizoens patronen te visualiseren
result = seasonal_decompose(df_time['cnt'], model='additive', period=24*7)
# Maken van de plot
plt.figure(figsize=(15, 5))
result.seasonal.plot()
plt.title('Seizoens Component')
plt.xlim('2012-04-01', '2012-05-01')
plt.xlabel('Datum')
plt.ylabel('Seizoens-waarde')
plt.show()
Wat hier te zien valt is dat er een duidelijke seizoensmatige beweging in de data zit. Momenteel is er een patroon dat zich 4 maal volledig toont. Vanuit deze wetenschap is er dan vanuit te gaan dat deze seizoenspatronen ongeveer een week duren. Om dit beter te kunnen bekijken zoomen we nu in op een periode van 10 dagen. Hierin zou het patroon zich moeten uiten en zal er een nieuwe herhaling aan het einde moeten beginnen.
# Gebruiken van seasonal_decompose om seizoens patronen te visualiseren
result = seasonal_decompose(df_time['cnt'], model='additive', period=24*7)
# Maken van de plot
plt.figure(figsize=(15, 5))
result.seasonal.plot()
plt.title('Seizoens Component')
plt.xlim('2012-04-01', '2012-04-10')
plt.xlabel('Datum')
plt.ylabel('Seizoens-waarde')
plt.show()
Dagelijks herhaalt zich een soortgelijk patroon. Er zijn kleine verschillen in aan te wijzen, maar de meeste dagen bevatten 2 pieken. Deze pieken vormen zich langzamerhand samen, waarna het weer scheidt naar twee pieken. De periode die het duurt om same te voegen is 5 dagen, terwijl de periode om te splitsen 2 dagen duurt. Het is aannemelijk dat het hierbij dus om de werkdagen en het weekend zal gaan.
# Berekenen aantal datapunten
N = len(df_time)
# Uitvoeren Fourier Transform
ft = np.fft.fft(df_time['cnt'])
# Bereken de kracht van de frequenties
magnitude = 2.0/N * np.abs(ft[:N//2])
# Zoeken en tonen van de pieken
peaks, _ = find_peaks(magnitude, height=75)
print(f"Peaks found at {peaks}")
# Plotten
plt.figure(figsize=(15,5))
plt.plot(magnitude)
plt.xlabel('Frequentie (1/N)')
plt.ylabel('Amplitude')
plt.title('Frequentie domein van Verhuursaantallen')
plt.plot(peaks, magnitude[peaks], ".", markersize=10)
plt.show()
Peaks found at [692 694]
Nu de locatie van de hoogste piek is gevonden, kunnen de frequentie en de periode worden berekend. Doordat peaks een lijst is met 2 waarden, kan er voor beide waarden de periode worden berekend.
freq = peaks * (1 / N)
periode = 1 / freq
print(f'De periode van FFT is: {np.round(periode, 0)}')
De periode van FFT is: [24. 24.]
Een periode geeft aan dat de Fourier Transform een periodieke herhaling ziet om de 24 tijdsperioden. In onze dataset betekent dit dat er elke dag een patroon plaats vind. Hoewel we al weten dat dit patroon kan verschillen op basis van weekend of geen weekend, was er ook te zien dat van maandag tot en met vrijdag de bewegingen van de grafiek gelijk aan elkaar waren. Om deze reden zullen de kolommen 'day_of_week' en 'is_weekend' worden togevoegd, zodat de patronen gebruikt kunnen worden voor Machine Learning.
Nu zal er worden gekeken naar de fourier features die toegevoegd kunnen worden.
# Aanmaken fourier voor W en D
fourier_w = CalendarFourier(freq="W", order=1)
fourier_d = CalendarFourier(freq="D", order=1)
dp = DeterministicProcess(
index=df_time.index,
constant=False,
order=1,
seasonal=False,
additional_terms=[fourier_w, fourier_d],
drop=True,
)
fourier_feat = dp.in_sample()
ax = fourier_feat.plot(y='sin(1,freq=D)',
xlim=['2012-01-01', '2012-01-03'])
ax = fourier_feat.plot(ax=ax, y='cos(1,freq=D)',
xlim=['2012-01-01', '2012-01-03'])
plt.show()
ax = fourier_feat.plot(y='sin(1,freq=W-SUN)',
xlim=('2012-01-01', '2012-01-15'))
ax = fourier_feat.plot(ax=ax, y='cos(1,freq=W-SUN)',
xlim=('2012-01-01', '2012-01-15'))
plt.show()
Zoals er te zien is in beide van de bovenstaande grafieken is er een mooie sinus beweging aanwezig. Deze sinusbeweging kan vervolgens als feature toegevoegd worden aan de datasets, omdat deze het verloop van de week of dag aangeeft.
Timeseries data is data verzameld over een periode van tijd. Deze data bestaat uit meerdere componenten, namelijk: de trendcomponent, seizoenscomponent, cyclische component en een willekeurig component. De trend component gaat over het gedrag van de data op de lange termijn, en negeert schommelingen op de korte termijn. Door het analyseren van de trendcomponent kunnen we data beter begrijpen, wat belangrijk is voor het voorspellen van toekomstige data. (Learn20stat, 2023)
Eerst onderzoeken we in hoeverre de trend component aanwezig is in onze data. Hier voor zijn er meerdere methoden die we kunnnen gebruiken. Een veel gebruikte methode is om moving averages te gebruiken.
df_trend = data_train.set_index('date_hour')
df_trend
| holiday | weathersit | temp | atemp | hum | windspeed | cnt | |
|---|---|---|---|---|---|---|---|
| date_hour | |||||||
| 2011-01-01 00:00:00 | 0 | 1 | 0.24 | 0.2879 | 0.81 | 0.000000 | 16 |
| 2011-01-01 01:00:00 | 0 | 1 | 0.22 | 0.2727 | 0.80 | 0.000000 | 40 |
| 2011-01-01 02:00:00 | 0 | 1 | 0.22 | 0.2727 | 0.80 | 0.000000 | 32 |
| 2011-01-01 03:00:00 | 0 | 1 | 0.24 | 0.2879 | 0.75 | 0.000000 | 13 |
| 2011-01-01 04:00:00 | 0 | 1 | 0.24 | 0.2879 | 0.75 | 0.000000 | 1 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 2012-11-30 19:00:00 | 0 | 1 | 0.32 | 0.3485 | 0.66 | 0.000000 | 377 |
| 2012-11-30 20:00:00 | 0 | 1 | 0.32 | 0.3485 | 0.66 | 0.000000 | 245 |
| 2012-11-30 21:00:00 | 0 | 1 | 0.30 | 0.3182 | 0.75 | 0.105325 | 183 |
| 2012-11-30 22:00:00 | 0 | 1 | 0.30 | 0.3333 | 0.75 | 0.000000 | 163 |
| 2012-11-30 23:00:00 | 0 | 2 | 0.30 | 0.3182 | 0.75 | 0.105325 | 110 |
16637 rows × 7 columns
# Index naar datetime
df_trend.index = pd.to_datetime(df_trend.index, format='%Y-%m-%d %H:%M:%S')
df_trend = df_trend.asfreq('H')
# Linear regression om de trendlijn te maken
X = pd.to_numeric(df_trend.index).values.reshape(-1, 1)
y = df_trend['cnt'].values
trendmodel = LinearRegression().fit(X, y)
# trendlijn voorspellen
trendline = trendmodel.predict(X)
# Originele timeseries plotten met trendlijn
plt.figure(figsize=(12, 6))
plt.plot(df_trend['cnt'], label='Original timeseries')
plt.plot(df_trend.index, trendline, label='Trendline', color='red')
plt.legend()
plt.title('Timeseries and Trendline')
plt.show()
In de grafiek is te zien dat er per jaar een boog gaat. De verhuuraantallen zijn in de zomer hoger. Goed te zien is dat de lijn in het tweede jaar in veel plekken hoger is dan in het vorige jaar. Dit kan een sterke opwaartse trend aangeven. Omdat de data over twee jaar gaat is het niet met zekerheid te zeggen of deze trend volhoud. bronnen voor dit hoofdstuk
Learn20stat. (2023, March 30). Components of time series. Learn Statistics. https://www.learn-stat.com/what-are-components-of-time-series/
Autocorrelatie is de correlatie/relatie van een tijdsreeks met een eerdere versie in tijd van hetzelfde dataframe. Het helpt om herhalende patronen te vinden.
Cycli verwijzen naar herhalende patronen of fluctuaties in de tijdreeksgegevens.
(Howell, 2023), (Learn Time Series Tutorials | Kaggle, z.d.)
# These are the "time lags"
shifts = np.arange(1, 25).astype(int)
# Use a dictionary comprehension to create name: value pairs, one pair per shift
shifted_data = {"lag_{}_hour".format(hour_shift): data_train['cnt'].shift(hour_shift) for hour_shift in shifts}
# Convert into a DataFrame for subsequent use
data_train_shifted_data = pd.DataFrame(shifted_data)
# Drop rows with NaN values introduced by the shift operation
data_train_shifted_data.dropna(inplace=True)
# Assuming 'target_variable' is the column you want to predict
target_variable = data_train['cnt'].iloc[data_train_shifted_data.index]
model = LinearRegression()
model.fit(data_train_shifted_data, target_variable)
fig, ax = plt.subplots()
ax.bar(data_train_shifted_data.columns, model.coef_)
ax.set(xlabel='Coefficient name', ylabel='Coefficient value')
# Set formatting so it looks nice
plt.setp(ax.get_xticklabels(), rotation=45, horizontalalignment='right')
# Plot the first 100 samples of each
ax = data_train_shifted_data.iloc[:100].plot(cmap=plt.cm.viridis)
target_variable.iloc[:100].plot(color='r', lw=2)
ax.legend(loc='best')
plt.show()
Een lag is het tijdsverschil tussen twee punten in een tijdreeks, en de autocorrelatie meet de mate van correlatie tussen deze twee punten. We hebben hier voor 24 lags gekozen want een dag bestaat uit 24 uur. In de barplot zie je de coefficiente waarde. en in de grafiek daar onder zie je de oorspronkelijke en voorspelde waarden voor de eerste 100 samples.
fig, (ax1, ax2) = plt.subplots(2,1, figsize=(12,8))
# Plot ACF
plot_acf(data_train['cnt'],lags=24, zero=False, ax=ax1)
ax1.set_ylim(0.2, -0.2)
# Plot PACF
plot_pacf(data_train['cnt'],lags=24, zero=False, ax=ax2)
ax2.set_ylim(0.2, -0.2)
plt.show()
Hier zie je autocorrelation function (ACF) en een partial autocorrelation function (PACF). De ACF meet de totale correlatie op een bepaalde vertraging en de PACF meet alleen de directe correlatie op een bepaalde vertraging. Uit de ACF kun je de Moving Average halen en uit de PACF kun je de AutoRegresive component halen. Deze grafieken gebruik onder andere wanneer je een (S)AR(I)MA(X) tijdreeks maakt. Hier verder in het bestand geven we meer uitleg over de ACF en PACF bij het onderdeel H.3.8SARIMAX. (Iamleonie, 2022)
result = seasonal_decompose(df_time['cnt'], model='additive', period=24*7)
# Maken van de plot
plt.figure(figsize=(15, 5))
result.seasonal.plot()
plt.title('Seizoens Component')
plt.xlim('2012-04-01', '2012-05-01')
plt.xlabel('Datum')
plt.ylabel('Seizoens-waarde')
plt.show()
Zoals hier boven in het kopje H2.3.1 Seizoens Analyse (met FFT) is te zien dat er een seizoens patroon is. Dit slaat ook op cycli want er zit een herhalend patroon in de tijdsreeks. Kijk maar naar de grafiek hierboven. Elke week herhaalt het zich weer.
Eerst zullen er de benodigde functies worden aangemaakt, om op die manier gemakkelijk de verschillende nieuwe time-series kolommen aan te maken. In de loop van het project zullen er kolommen worden toegevoegd aan deze functie.
Dit zijn de kolommen waarmee we zijn begonnen:
Hier is een korte omschrijving over de later toegevoegde kolommen en waarom ze zijn toegevoegd:
def data_voorbereiding(base_df, dt_kol):
"""
Een functie om een csv in te lezen en verschillende
toepassingen uit te voeren op aanwezige timeseries
data.
Parameters:
----------
base_df : pd.DataFrame
De naam van het dataframe waar de extra
timeseries kolommen worden gemaakt.
dt_kol : str
De naam van de kolom met de Timeseries data.
Returns:
----------
df : pd.DataFrame
Een pandas.DataFrame waarin verschillende extra
kolommen zijn toegevoegd waarin timeseries is
toegepast.
"""
# Omzetten van date_hour naar datetime
df = base_df.copy()
df[f'{dt_kol}'] = pd.to_datetime(df[f'{dt_kol}'])
# Toegevoegd op basis van eigen ingeving
df['Jaar'] = df[f'{dt_kol}'].dt.year.astype(int)
df['Maand'] = df[f'{dt_kol}'].dt.month.astype(int)
df['Uur'] = df[f'{dt_kol}'].dt.hour.astype(int)
# Toegevoegd na aanleiding FFT
df['is_weekend'] = (df[f'{dt_kol}'].dt.dayofweek >= 5).astype(int)
# Aanmaken fourier voor W en D
fourier_w = CalendarFourier(freq="W", order=1)
fourier_d = CalendarFourier(freq="D", order=1)
# Index datetime maken
base_df[f'{dt_kol}'] = pd.to_datetime(base_df[f'{dt_kol}'])
base_df = base_df.set_index(f'{dt_kol}')
# Uitvoeren van DeterministicProcess
dp = DeterministicProcess(
index=base_df.index,
constant=False,
order=1,
seasonal=False,
additional_terms=[fourier_w, fourier_d],
drop=True,
)
# Aanmaken van fourier dataframe
fourier_feat = dp.in_sample()
# Mergen van dataframes voor alle data
df_fou = pd.merge(
df, fourier_feat[['sin(1,freq=D)', 'cos(1,freq=D)',
'sin(1,freq=W-SUN)', 'cos(1,freq=W-SUN)']],
left_on='date_hour', right_index=True)
return df_fou
Deze functies zullen helpen om nieuwe kolommen toe te voegen aan data_train, en later ook aan data_test.
# Toepassen data_voorbereiding
df = data_voorbereiding(data_train, 'date_hour')
df.head()
| date_hour | holiday | weathersit | temp | atemp | hum | windspeed | cnt | Jaar | Maand | Uur | is_weekend | sin(1,freq=D) | cos(1,freq=D) | sin(1,freq=W-SUN) | cos(1,freq=W-SUN) | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2011-01-01 00:00:00 | 0 | 1 | 0.24 | 0.2879 | 0.81 | 0.0 | 16 | 2011 | 1 | 0 | 1 | 0.000000 | 1.000000 | -0.974928 | -0.222521 |
| 1 | 2011-01-01 01:00:00 | 0 | 1 | 0.22 | 0.2727 | 0.80 | 0.0 | 40 | 2011 | 1 | 1 | 1 | 0.258819 | 0.965926 | -0.982566 | -0.185912 |
| 2 | 2011-01-01 02:00:00 | 0 | 1 | 0.22 | 0.2727 | 0.80 | 0.0 | 32 | 2011 | 1 | 2 | 1 | 0.500000 | 0.866025 | -0.988831 | -0.149042 |
| 3 | 2011-01-01 03:00:00 | 0 | 1 | 0.24 | 0.2879 | 0.75 | 0.0 | 13 | 2011 | 1 | 3 | 1 | 0.707107 | 0.707107 | -0.993712 | -0.111964 |
| 4 | 2011-01-01 04:00:00 | 0 | 1 | 0.24 | 0.2879 | 0.75 | 0.0 | 1 | 2011 | 1 | 4 | 1 | 0.866025 | 0.500000 | -0.997204 | -0.074730 |
# Opnieuw aanmaken Grafieken class
graf = Grafieken(df, 'cnt')
# Toepassen timeplot method
graf.timeplot('Jaar', 'weathersit')
Zoals hierboven te zien is, zijn er in 2012 meer verhuursaantallen. Ook is er heel duidelijk te zien dat er een duidelijk verschil is tussen de verschillende weersituaties. Dat de weersomstandigheden een grote invloed hebben is mogelijk als het te verhuren product gebruikt wordt zonder afdekking van regen of andere neerslag. De boxplots tonen dat er geen outliers bij situatie 4 aanwezig zijn.
Bij een weerssituatie van klasse van 1 is er een algemeen hoger verhuur, dit neemt af naarmate je omhoog gaat op de schaal. Hieronder is verduidelijking over de schaal van deze kolom:
weathersit:
Voor verder onderzoek wordt er nu gekeken naar de verdeling van het verhuur per uur.
# Toepassen timeplot method
graf.timeplot('Uur', 'is_weekend')
De grafiek over de verhuuraantallen per uur geeft ook een duidelijke boodschap. Rond de uren van 6 t/m 8 vindt er een sterke stijging plaats, wat niet heel raar is aangezien de meeste mensen tussen 6 en 8 wakker worden. Rond de uren van 17 en 18 is er ook een gigantische piek. Echter is dit alleen van maandag t/m vrijdag, dit zijn de algemene werkdagen.
In de weekenden zit de hoogste piek rond het begin van de middag. Wat aan dit verschil ook opvalt is dat er in de daluren (behalve van 20 t/m 22) in het weekend altijd meer verhuur is. Terwijl in de spitsuren het juist drukker is op werkdagen.
Om iets meer inzicht te krijgen over het algemene verhuur gedrag is deze grafiek ook in de vorm van boxplots gemaakt, om zo de verdeling en de outliers ook duidelijk te zien. Het voordeel aan boxplots is dat er per waarde op de x-as een verdeling wordt getoond. Deze verdeling laten zien waar het gemiddelde ligt en of er outliers aanwezig zijn. Wat opvalt is de afwezigheid van outliers in de drukste uren van het verhuren en dat outliers ook enigzins afhankelijk zijn van tijd en de weekdag. Verder is het ook opvallend dat er geen outliers zijn aan de onderkant van de boxplots.
Tijdens deze stap gaat er gemodelleerd worden. Bij elk model zijn er ook tussenstappen waarbij features worden toegevoegd of verwijderd, om op deze manier de invloed van de features te kunnen beoordelen.
df_test = data_voorbereiding(data_test, 'date_hour')
df_test = df_test.set_index('date_hour').asfreq('H')
df_test
| holiday | weathersit | temp | atemp | hum | windspeed | Jaar | Maand | Uur | is_weekend | sin(1,freq=D) | cos(1,freq=D) | sin(1,freq=W-SUN) | cos(1,freq=W-SUN) | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| date_hour | ||||||||||||||
| 2012-12-01 00:00:00 | 0 | 1 | 0.26 | 0.3030 | 0.81 | 0.0000 | 2012 | 12 | 0 | 1 | 0.000000 | 1.000000 | -0.974928 | -0.222521 |
| 2012-12-01 01:00:00 | 0 | 1 | 0.26 | 0.3030 | 0.81 | 0.0000 | 2012 | 12 | 1 | 1 | 0.258819 | 0.965926 | -0.982566 | -0.185912 |
| 2012-12-01 02:00:00 | 0 | 2 | 0.26 | 0.3030 | 0.81 | 0.0000 | 2012 | 12 | 2 | 1 | 0.500000 | 0.866025 | -0.988831 | -0.149042 |
| 2012-12-01 03:00:00 | 0 | 2 | 0.26 | 0.2727 | 0.81 | 0.1343 | 2012 | 12 | 3 | 1 | 0.707107 | 0.707107 | -0.993712 | -0.111964 |
| 2012-12-01 04:00:00 | 0 | 1 | 0.26 | 0.2879 | 0.81 | 0.0896 | 2012 | 12 | 4 | 1 | 0.866025 | 0.500000 | -0.997204 | -0.074730 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2012-12-19 19:00:00 | 0 | 1 | 0.38 | 0.3939 | 0.50 | 0.3881 | 2012 | 12 | 19 | 0 | -0.965926 | 0.258819 | 0.593820 | -0.804598 |
| 2012-12-19 20:00:00 | 0 | 1 | 0.36 | 0.3485 | 0.57 | 0.2239 | 2012 | 12 | 20 | 0 | -0.866025 | 0.500000 | 0.563320 | -0.826239 |
| 2012-12-19 21:00:00 | 0 | 1 | 0.34 | 0.3182 | 0.61 | 0.2239 | 2012 | 12 | 21 | 0 | -0.707107 | 0.707107 | 0.532032 | -0.846724 |
| 2012-12-19 22:00:00 | 0 | 1 | 0.34 | 0.3485 | 0.61 | 0.0896 | 2012 | 12 | 22 | 0 | -0.500000 | 0.866025 | 0.500000 | -0.866025 |
| 2012-12-19 23:00:00 | 0 | 1 | 0.32 | 0.3333 | 0.66 | 0.1343 | 2012 | 12 | 23 | 0 | -0.258819 | 0.965926 | 0.467269 | -0.884115 |
456 rows × 14 columns
def submissie(model, suffix):
"""
Deze functie maakt de voorspelling op de test dataset
en vormt deze om tot een csv bestand om in te kunnen
leveren op Kaggle.
Parameters:
----------
model : Machine Learning model
De naam die is gegeven aan het ML-model dat
wordt gebruikt om te voorspellen
suffix : str
De laatste 'tag' voor de naam van het csv
bestand, zodat deze makkelijk te identificeren
is na de submission.
Returns:
----------
None
In plaats van een return maakt het een bestand
aan in de map Kaggle Submissions. De namen zien
er als volgt uit:
"""
# Fitten, voorspellen en veranderen naar integer
y_pred = model.fit(X, y).predict(df_test)
y_pred = y_pred.astype(int)
# Aanmaken df met alleen p_id en Outcome
test_predictions_df = pd.DataFrame(
{'date_hour': data_test['date_hour'],
'cnt': y_pred})
# Aanmaken van tijd
tijd = datetime.datetime.now().strftime("%m%d%H%M%S")
# Aanmaken csv bestand met timestamp
test_predictions_df.to_csv(
f'Kaggle Submissions/vs_{suffix}_{tijd}.csv',
index=False)
# print voor conformatie
print(f'vs_{suffix}_{tijd}.csv has been saved!')
def model_score(model):
"""
Een functie die de scores voor een model
berekend en deze toont.
Parameters:
----------
model : elk ML model
De naam van het gemaakte ML model
Returns:
----------
f1 : int
De f1 score van het gegeven model
"""
# Fitten en voorspellen met model
y_pred = model.fit(X_train, y_train).predict(X_test)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
# Tonen van de resultaten
return print(f'Root Mean Squared Error: {rmse}')
def grid_score(estimator, param_grid):
"""
Een functie die de scores voor een grid search
berekend en deze toont. Hierbij worden ook de
optimale parameters getoont.
Parameters:
----------
estimator : ML model
Het model waarop de grid search wordt toegepast
param_grid : dict
Een dictionary waarbij de keys parameter namen
zijn van het model en de values verschillende
waarden die de parameter aan kan nemen.
Returns:
----------
None
Returned twee prints met de optimale parameters
en beste f1 score
"""
# Uitvoeren van de grid search
gs = GridSearchCV(estimator=estimator,
param_grid=param_grid,
cv=cv,
scoring='neg_root_mean_squared_error',
n_jobs=-1)
# Fitten van de grid search
gs.fit(X_train, y_train)
# Tonen van de beste score en parameters
print(f"Beste RMSE: {abs(gs.best_score_)}")
print(f"Beste parameters:\n{gs.best_params_}")
def feature_scan(model):
"""
Een functie die een grafie met de feature importances kan maken.
Dit werkt in ieder geval voor de RF en XGB modellen.
Parameters:
----------
model : ML-model
Het model waarvan je de feature importances wilt tonen.
Returns:
----------
None
Laat een grafiek zien met daarin de invloed van alle features.
"""
# Aanmaken en sorteren invloed kolommen
invloed = pd.Series(model.feature_importances_, index=X.columns)
invloed = invloed.sort_values(ascending=True)
# Plotten grafiek met invloed
invloed.plot(kind='barh', figsize=(10, 6))
plt.ylabel('Features')
plt.xlabel('Invloed')
plt.show()
# Aanmaken SEED voor random_state
SEED = 42
# Aanmaken van X en y
df_train = df.set_index('date_hour').asfreq('H')
X = df_train.drop('cnt', axis=1)
y = df_train['cnt']
# Splitten van de data in train en test
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=456, random_state=SEED)
# Aanmaken van de folds
cv = TimeSeriesSplit(n_splits=5, test_size=456)
De waarde van een afhankelijke variabele voorspellen op basis van een of meer onafhankelijke variabelen. De formule van lineare regressie is y = b0 + b1x. b0 is het snijpunt met de y-as, b1 is de helling van de lijn, x is de onafhankelijke variabele en y is de afhankelijke variabele. (pp uit de lessen)
Een loss functie is een manier om te meten hoe goed een lineair regressiemodel de gegevens beschrijft. Het is een functie die de fout of afwijking tussen de werkelijke uitvoer en de voorspelde uitvoer kwantificeert. Een veel gebruikte loss functie voor lineaire regressie is de mean squared error (MSE) functie, die het gemiddelde van de kwadraten van de fouten berekent. Hoe kleiner de MSE, hoe beter het model de gegevens past. Voor deze opdracht gebruiken wij de root mean squared error (RMSE).
De formule is als volgt: $RMSE = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (y_i - \hat{y_i})^2}$
waar N: aantal waarnemingen is, i: de index, y_i: true labels en ^y_i: predicted labels is
(notebook loss_functions)
De vorm van regularisatie die van toepassing is op lineaire regression is de L1 en L2 regularisatie. Ook wel de Lasso (L1) en ridge (L2) regularisatie genoemd. Door regularisatie toe te passen, wordt een complex model tijdens het trainen vereenvoudigd. L1 voegt de “absolute waarde van de grootte” (“absolute value of magnitude”) van de coëfficiënt als penalty term toe aan de loss-functie. L2 voegt de “kwadratische omvang” (“squared magnitude”) van de coëfficiënt als penalty term toe aan de loss-functie. (Nagpal, 2022), (opdracht 1 diabetes)
Lasso (L1)
Wiskundige formule voor L1 is:
$\lambda||w||_1$
Waarbij:
Lasso (ook bekend als L1-regularisatie) voegt een term toe aan de loss functie die in gelijke verhouding staat tot de absolute waarden van de gewichten. Het effect dat Lasso heeft is dat het sommige gewichten reduceert to nul. Hierdoor worden sommige kenmerken volledig genegeeerd, waardoor feature selectie mogelijk wordt gemaakt. Door middel van feature selectie wordt de kans op overfitting verminderd. (Jain, 2023; Team, 2020), (opdracht 1 diabetes)
Ridge (L2)
De Wiskundige formule voor L2 is:
$\lambda w^2$
Waarbij:
Ridge (ook bekend als L2-regularisatie) voegt een term toe aan de loss functie die in gelijke verhouding staat tot de kwadraten van de gewichten. Het effect dat Ridge heeft is dat het de waarden van gewichten dicht naar de nul brengt, maar niet nul maakt. De ridge manier gaat de mogelijke overfitting tegen door te zorgen dat de gewichten geen te grootte waarden kunnen aannemen. (Jain, 2023; Team, 2020), (opdracht 1 diabetes)
Nu het duidelijk is wat lineaire Regressie inhoud, kan er een model worden opgesteld.
# Opstellen van Linear Regression model
lr = LinearRegression(n_jobs=-1)
# Bekijken van de RMSE
model_score(lr)
Root Mean Squared Error: 117.2269334691733
Nu zal er een parameter grid aangemaakt worden voor een Grid Search CV, om zo de optimale hyperparameters te kunnen bepalen.
param_lr = {
'fit_intercept': [True, False],
'copy_X' : [True, False],
'positive' : [False, True]
}
# Uitvoeren van Grid Search CV
grid_score(lr, param_lr)
Beste RMSE: 126.75350985874736
Beste parameters:
{'copy_X': True, 'fit_intercept': True, 'positive': False}
Deze parameters kunnen nu worden ingevoerd in het model.
# Toepassen van optimale parameters
lr_tuned = LinearRegression(n_jobs=-1)
# Bekijken van de RMSE
model_score(lr_tuned)
Root Mean Squared Error: 117.2269334691733
Nu het model is geoptimaliseerd wordt er met de submissie functie een inleverbaar bestand gemaakt voor op Kaggle.
#submissie(lr_tuned, 'LR')
Wij hebben gekozen om een k-nearest neighbors model te gebruiken. Een K-nearest neighbours model (KNN) kan gebruikt worden voor classificatie of regressie problemen. In een KNN model word voor het voorspellen gekeken naar datapunten die in de buurt van de waarde die we willen voorspellen zitten, ofwel de "neighbours". Het model neemt dus de X waarden van de data die voorspeld moet worden, en vergelijkt dit met de X waarden van de trainingsdata die er vergelijkbaar uit ziet. Zo kan er een voorspelling worden gemaakt over de y waarde. Hier voor is er de mogelijkheid gebruik te maken van 3 verschillende metrics.(IBM, 2023)
Er zijn enkele voordelen verbonden aan een KNN model. Zo is het model niet zo complex als andere modellen, wat de implementatie makkelijker maakt. Verder zijn er weinig hyperparameters, dit is zowel een voordeel als een nadeel.
Omdat het model redelijk simpel is, wordt er op andere vlakken ingeleverd. Zo kan het model met grotere datasets veel eisen van de computer. Ook is er een hogere kans op overfitting. (GeeksforGeeks, 2023)
# opstellen K-Nearest Neighbors model
knn = KNeighborsRegressor(n_jobs=-1)
model_score(knn)
Root Mean Squared Error: 53.38359939354527
Aan de hand van de onderstaande code kan er door trial en error een optimale waarde voor het aantal neighbours worden bepaald.
# Bekijken aantal neighbors
train_accuracies = {}
test_accuracies = {}
neighbors = range(1, 11)
for neighbor in neighbors:
knn = KNeighborsRegressor(n_neighbors=neighbor)
knn.fit(X_train, y_train)
train_accuracies[neighbor] = knn.score(X_train, y_train)
test_accuracies[neighbor] = knn.score(X_test, y_test)
display(train_accuracies, test_accuracies)
{1: 0.9999815166962486,
2: 0.9643702994604253,
3: 0.95209467887222,
4: 0.947853355542497,
5: 0.9451850853115388,
6: 0.9418743474241642,
7: 0.9390014376078347,
8: 0.9362817072430636,
9: 0.9333881860709422,
10: 0.9310182177427528}
{1: 0.8355177220420447,
2: 0.8832341262210887,
3: 0.8887120763203489,
4: 0.8936747310379309,
5: 0.897679396895329,
6: 0.9028353960226173,
7: 0.9100914424354014,
8: 0.9107787975884452,
9: 0.907995193454967,
10: 0.9083116354148302}
Door het optimale aantal neighbors in te voeren kan het model verbeterd worden.
knn = KNeighborsRegressor(n_neighbors=4)
model_score(knn)
Root Mean Squared Error: 54.4182475724793
Nu zal er een parameter grid aangemaakt worden voor een Grid Search CV, om zo de optimale hyperparameters te kunnen bepalen.
# parameters voor de GridSearch
params_knn = {'weights':['uniform', 'distance'],
'metric':['minkowski','euclidean','manhattan' ],
'algorithm':['ball_tree', 'kd_tree', 'brute']}
grid_score(knn, params_knn)
Beste RMSE: 55.75803076569273
Beste parameters:
{'algorithm': 'ball_tree', 'metric': 'minkowski', 'weights': 'uniform'}
Deze parameters kunnen nu worden ingevoerd in het model.
# uitkomsten GridSearch toepassen
knn_tuned = KNeighborsRegressor(algorithm='ball_tree',
n_jobs=-1,
metric='minkowski',
n_neighbors=4,
weights='uniform')
model_score(knn_tuned)
Root Mean Squared Error: 54.4182475724793
De Decision Tree Regressor wordt, in tegenstelling tot de Decision Tree Classifier, gebruikt voor regressieproblemen. Een Decision Tree is opgebouwd uit verschillende knooppunten (nodes). De wortel (root node) is het startknooppunt van de boom en bevat alle beschikbare data. De wortel wordt vervolgens gesplitst in twee interne knooppunten (internal nodes) op basis van een bepaald kenmerk en drempelwaarde. De boom wordt gesplitst op basis van een bepaald kenmerk en drempelwaarde. De keuze van welk kenmerk en welke drempelwaarde te gebruiken wordt bepaald door criteria zoals de Gini-index of de entropie. Deze criteria meten de homogeniteit van de gegevens aan elke kant van de splitsing. Elke interne knoop vertegenwoordigt een beslissingspunt waar de boom bepaalt welke richting te volgen op basis van de kenmerken van de gegeven datapunten. Elk intern knooppunt wordt op zijn beurt weer gesplitst in twee, hetzij in andere interne knooppunten, hetzij in bladeren (leaf nodes). Deze splitsingen worden herhaald totdat de boom de bladeren bereikt. Een blad is het uiteindelijke niveau van de boom en vertegenwoordigt de voorspelde uitkomst voor de gegeven invoer.
(Classification and regression trees - Decision tree for classification, z.d.), (K, 2021)
# Opstellen van Dicision Tree Regressor model
dt = DecisionTreeRegressor(random_state=42)
# Bekijken van de RMSE
model_score(dt)
Root Mean Squared Error: 54.01557199126803
Nu zal er een parameter grid aangemaakt worden voor een Grid Search CV, om zo de optimale hyperparameters te kunnen bepalen.
param_dtr = {
'criterion': ['squared_error', 'friedman_mse'],
'max_depth': [None, 10, 20, 30],
'min_samples_split': [2, 5],
'min_samples_leaf': [1, 2, 4],
}
# Uitvoeren van Grid Search CV
grid_score(dt, param_dtr)
Beste RMSE: 56.58902840386772
Beste parameters:
{'criterion': 'squared_error', 'max_depth': 20, 'min_samples_leaf': 4, 'min_samples_split': 2}
Met deze parameters kan Decision Tree opnieuw worden uitgevoerd. Dit zal een beter resultaat leveren dan er voorheen was te zien.
# Toepassen van optimale parameters
dt_tuned = DecisionTreeRegressor(criterion='squared_error',
max_depth=None,
min_samples_leaf=4,
min_samples_split=2)
# Bekijken van de RMSE
model_score(dt_tuned)
Root Mean Squared Error: 49.20313045211243
Nu het getunede model van Decision Tree is uitgevoerd, kan er gekeken worden naar de verschillende feature importances. Een grafiek met deze data geeft aan welke kolommen invloed hebben en hoe sterk de invloed is.
feature_scan(dt_tuned)
Nu het model is geoptimaliseerd wordt er met de submissie functie een inleverbaar bestand gemaakt voor op Kaggle.
#submissie(dt_tuned, 'DT')
Random forest begint met het aanmaken van meerdere decision tree modellen, die elk tegelijk gaan trainen en voorspellen. Deze methode heet ensemble learning, het gebruiken van verschillende of meerdere modellen om tot een voorspelling te komen. Bij een Random Forest wordt er enkel gebruik gemaakt van Decision Trees. Nadat elk van de trees een voorspelling heeft gemaakt, wordt uit die voorspellingen het gemiddelde berekent om zo tot een antwoord te komen op het regressie probleem.
Het aantal gebruikte features is bij een Random Forest afhankelijk van een hyper-parameter: max_features. Deze parameter zorgt ervoor dat er een limiet zit op de features, zodat het ensemble model niet te afhankelijk kan worden van een paar features. Ook niet alle data wordt bij elke tree ingezet, om te zorgen dat er een extra willekeur inzit om overfitting tegen te gaan. (Chakure, 2022)
Deze eigenschappen van Random Forest zijn waarom wij het model gaan gebruiken.
# Opstellen van Random Forest Regressor model
rf = RandomForestRegressor(random_state=SEED,
n_jobs=-1)
# Bekijken van de RMSE
model_score(rf)
Root Mean Squared Error: 37.84325840402451
Nu zal er een parameter grid aangemaakt worden voor een Grid Search CV, om zo de optimale hyperparameters te kunnen bepalen.
param_rf = {
'n_estimators' : [50, 100, 200],
'max_depth' : [10, 20, 30],
'min_samples_leaf' : [1, 2, 4]
}
grid_score(rf, param_rf)
Beste RMSE: 45.310138190556685
Beste parameters:
{'max_depth': 30, 'min_samples_leaf': 1, 'n_estimators': 200}
Nu deze parameters bekend zijn, kunnen deze worden ingevoerd in het model. Op deze manier zal het model beter presteren.
# Toepassen van optimale parameters
rf_tuned = RandomForestRegressor(n_estimators=200,
max_depth=30,
min_samples_leaf=1,
random_state=SEED,
n_jobs=-1)
# Bekijken van de RMSE
model_score(rf_tuned)
Root Mean Squared Error: 37.78426103437784
Nu het getunede model van Random Forest is uitgevoerd, kan er gekeken worden naar de verschillende feature importances. Een grafiek met deze data geeft aan welke kolommen invloed hebben en hoe sterk de invloed is.
feature_scan(rf_tuned)
De bovenstaande grafiek toont dat, voor het maken van een random forest voorspelling, de waarden van uur en temperatuur het belangrijkste zijn. De kolom holiday heeft duidelijk het minste invloed op de voorspelling. Verder is uur duidelijk de waarde waarop het model haar voorspellingen het meest heeft gebaseerd. In de EDA was er ook duidelijk te zien dat het uur van de dag grote invloed had op de verhuuraantallen, dus dat het model hier voor heeft gekozen is logisch.
#submissie(rf_tuned, 'RF_2fa')
Dit model heeft enkele voordelen die er toe hebben geleid om dit model te kiezen. Volgens Kumar (2019) is een van de voordelen de regularisatie van een XGBoost model. Dit model bevat zowel Lasso (L1) als Ridge (L2) regularisatie, wat ervoor zorgt dat het model moeilijker overfit. Via de Scikit Learn libary is het mogelijk om verschillende waarden, als hyperparameters, toe te voegen aan het model. De parameter alpha is voor L1 regularisatie en de parameter lambda is voor L2 regularisatie.
Een ander voordeel is de manier waarop het model pruning toepast. Bij een XGBoost model worden de bomen pas aan het einde van het model gepruned. Dit zorgt ervoor dat als er een enkele split met een negative loss voorkomt, het model niet de tak daar laat stoppen. Het model zal eerst doorgaan tot de aangegeven max_depth voordat deze gaat prunen. (Kumar, 2019)
Daarnaast is een ander groot voordeel dat XGBoost een parallel lopende ensemble is. Dit betekent dat de bomen tegelijkertijd kunnen worden opgebouwd, in tegen stelling to achtereenvolgend. Dit zorgt voornamelijk in een kortere run time bij het trainen van het model. (Hanchman, 2023)
Als eerste stap trainen we een model met enkel de standaard parameters.
# Opstellen van XGBoost Regressor model
xgb = XGBRegressor(n_jobs=-1,
random_state=SEED)
# Bekijken van de RMSE
model_score(xgb)
Root Mean Squared Error: 37.458258453583944
Nu zal er een parameter grid aangemaakt worden voor een Grid Search CV, om zo de optimale hyperparameters te kunnen bepalen.
# Opstellen van de parameter dict
param_xgb = {
'n_estimators': [100, 200, 300],
'learning_rate': [0.01, 0.1, 1.0],
'max_depth': [5, 10, 20]
}
# Uitvoeren van de gridsearch
grid_score(xgb, param_xgb)
Beste RMSE: 41.866621142492704
Beste parameters:
{'learning_rate': 0.1, 'max_depth': 5, 'n_estimators': 300}
Nu deze parameters bekend zijn, kunnen deze worden ingevoerd in het model. Op deze manier zal het model beter presteren.
# Toepassen van optimale parameters
xgb_tuned = XGBRegressor(n_estimators=300,
learning_rate=0.1,
max_depth=5,
random_state=SEED,
n_jobs=-1)
# Bekijken van de RMSE
model_score(xgb_tuned)
Root Mean Squared Error: 36.402051993205
Ook bij XGBoosting zal worden gekeken naar de feature importance.
feature_scan(xgb_tuned)
De bovenstaande grafiek toont dat, voor het maken van een XGBoost voorspelling, de waarden van uur en jaar het belangrijkste zijn. De kolom windspeed heeft het minste invloed op de voorspelling. De waarden van uur en jaar liggen dichtbij elkaar qua invloed. Wat bij deze grafiek ook duidelijk is te zien, is dat alle waarden iets hoger liggen dan bij random forest. Nu kan er een bestand worden gemaakt om op Kaggle in te leveren.
# submissie(xgb_tuned, 'XGB_new')
Binnen machine learning zijn er vele verschillende ensemble modellen, dit zijn modellen die gebruik maken van andere modellen om zo de voorspellingte kunnen vormen. Binnen het kader van Ensemble learning zijn er drie verschillende typen: bagging, boosting en stacking.
Bagging:
Bij bagging wordt de data gesplitst in meerdere gelijke delen, een voorbeeld hiervan komt voor bij RandomForestRegressor. Elke tree die RandomForest aanmaakt gebruikt een andere subset aan features om de voorspelling te maken. Van deze voorspelling wordt vervolgens een gemiddelde berekent, deze waarde dient als de voorspelling van het model.
Boosting:
Bij boosting wordt er eerste een model getrained op de dataset. Vervolgens gaat het tweede model zich bezig houden met het verbeteren van de voorspelling, dit gebeurt door de gewichten op de verkeerd voorspelde datapunten te vergroten. Dit herhaalt zich totdat er geen verbetering meer plaats vind of tot het maximale aantal rondes is bereikt. Een populair voorbeeld van boosting in ML is het gebruik van ADABoost.
Stacking:
Als keuze voor een ensemble model is er gekozen om te werken met StackingRegressor. Net als bij bagging en boosting combineerd stacking meerdere modellen om tot een voorspelling te komen voor de dataset. Anders dan bij bagging zijn er verschillende modellen die op dezelfde data worden getrained. En anders dan bij boosting word een ander model gebruikt om deze voorspellingen om te zetten naar een enkele waarde (Brownlee, 2021). Dit zorgt ervoor dat een stacking model de voor en nadelen van meerdere modellen kan combineren om een, hopelijk, sterker model te bouwen.
Wiskundig gezien kan je elke collectie van voorspelling als volgt classificeren: $\hat{y}_i,_s$. Dit staat voor de voorspelling $\hat{y}$ gebaseerd op de basis modellen $i$ op basis van sample $s$. Deze waarden worden als feature gebruikt voor het trainen van het zogenoemde meta-model, het model dat de uiteindelijke voorspellingen maakt. Bij het voorspellen van de target word er nieuwe data door de basis modellen gevoerd. Deze modellen geven een voorspelling, zodra deze voorspellingen zijn gemaakt gaat het meta-model voorspellen.
# Opzetten basis- en meta-modellen
basis_modellen = [
('knn', knn_tuned),
('rf', rf_tuned),
('lr', lr_tuned)
]
meta_model = xgb_tuned
# Toepassen van StackingRegressor
sr = StackingRegressor(estimators=basis_modellen,
final_estimator=meta_model,
n_jobs=-1)
# Bekijken van de RMSE
model_score(sr)
Root Mean Squared Error: 40.11728088527348
Nu zal er een parameter grid aangemaakt worden voor een Grid Search CV, om zo de optimale hyperparameters te kunnen bepalen.
# Opstellen van parameter dict
param_sr = {
'passthrough' : [True, False],
'cv' : [5, 10, None]
}
# Uitvoeren van gridsearch
grid_score(sr, param_sr)
Beste RMSE: 42.15045747739301
Beste parameters:
{'cv': 5, 'passthrough': True}
Nu deze parameters bekend zijn, kunnen deze worden ingevoerd in het model. Op deze manier zal het model beter presteren.
# Opzetten basis- en meta-modellen
basis_modellen = [
('knn', knn_tuned),
('rf', rf_tuned),
('lr', lr_tuned)
]
meta_model = xgb_tuned
# Toepassen van StackingRegressor
sr_tuned = StackingRegressor(estimators=basis_modellen,
final_estimator=meta_model,
cv=5,
passthrough=True,
n_jobs=-1)
# Bekijken van de RMSE
model_score(sr_tuned)
Root Mean Squared Error: 37.7191700700007
Nu het model is geoptimaliseerd wordt er met de submissie functie een inleverbaar bestand gemaakt voor op Kaggle.
# submissie(sr_tuned, 'SR')
Hybride Machine Learning is het gebruiken van meerdere modellen om te zorgen voor een betere voorspelling. Vaak wordt er een combinatie gebruikt van een simpel (bijvoorbeeld Lineare Regressie) en een meer complex (bijvoorbeeld Random Forest) model (Holbrook, 2023). Tijdens het uitvoeren van de code worden de nodige stappen voor het model uitgelegd.
Om te beginnen met de code zal het eerste model, Linear Regression, de voorspellingen maken.
# Eerste model: Lineaire Regressie
lr_hm = LinearRegression(n_jobs=-1).fit(X_train, y_train)
# Voorspellen van de waarden voor train en test
y_pred_lr = lr_hm.predict(X_train)
pred_test = lr_hm.predict(X_test)
Nu de eerste voorspellingen zijn gemaakt, kan het tweede model aan de slag gaan. Bij een hybride model werkt dit anders dan bij een ensemble model. Waar er bij een ensemble gebruikt wordt gemaakt van het gemiddelde van de voorspellingen, zal een hybride model gebruik maken van de residuals van het eerste model.
# Berekenen van de residuals
residuals = y_train - y_pred_lr
Deze nieuwe variabele bevat het verschil tussen de werkelijke waarden en de voorspelde waarden. Het tweede model zal deze verschillen proberen te voorspellen, om deze vervolgens bij de andere voorspellingen op te kunnen tellen. Dit tweede model, XGBoost, zal dus aan de hand van de residuals gaan voorspellen wat de verschillen zullen zijn tussen de twee waarden.
# Tweede model: XGBoost Regressor
xgb_hm = XGBRegressor(n_estimators=300,
learning_rate=0.1,
max_depth=5,
random_state=SEED,
n_jobs=-1).fit(X_train, residuals)
# Voorspellen van de residuals
y_pred_xgb = xgb_hm.predict(X_test)
Nu de verschillen zijn voorspeld kunnen deze bij de voorspelde waarden worden gevoegd. Op deze manier komen we bij de daadwerkelijke voorspelling komen van het hybride model.
# Samenvoegen van voorspellingen
pred = np.round(pred_test + y_pred_xgb).astype(int)
# Evalueren van hybride model
mse = mean_squared_error(y_test, pred)
print(f'Root Mean Squared Error: {np.sqrt(mse)}')
Root Mean Squared Error: 36.66069528728783
Nu het laatste model is getrained kunnen de feature importances worden bekeken.
feature_scan(xgb_hm)
Nu de feature importances bekeken zijn en de RMSE is berekend kunnen de waarden omgezet worden tot een dataframe die ingeleverd kan worden op kaggle. Deze code is gebaseerd op de code uit de submissie functie, echter was het niet mogelijk om de functie te gebruiken. Dit is wegens de samenvoeging van de voorspellingen.
# Voorspellen op df_test
lr_test = lr_hm.predict(df_test)
xgb_test = xgb_hm.predict(df_test)
# Samenvoegen voorspelling
y_pred = np.round(lr_test + xgb_test).astype(int)
# Aanmaken df met alleen p_id en Outcome
test_predictions_df = pd.DataFrame(
{'date_hour': data_test['date_hour'],
'cnt': y_pred})
# Aanmaken van tijd
# tijd = datetime.datetime.now().strftime("%m%d%H%M%S")
# Aanmaken csv bestand met timestamp
# test_predictions_df.to_csv(
# f'Kaggle Submissions/vs_HML_try_{tijd}.csv',
# index=False)
# print voor conformatie
# print(f'vs_HML_try_{tijd}.csv has been saved!')
De informatie en code hier onder heb ik uit verschillende bronnen gehaald:
Sarimax is een tijdserie model en de afkorting staat voor Seasonal AutoRegresive Integrated Moving Avarege eXogenous.
Als eerst gaan we een Dickey-Fuller test uitvoeren. Deze bepaalt of een dataset stationair is of niet. De data is stationair als de waarde van p kleiner of gelijk is aan het significantie level. Als de p-waarde groter is, is de data niet stationair
#Augmented Dickey–Fuller test:
print('Results of Dickey Fuller Test:')
dftest = adfuller(df_train['cnt'], autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
print(dfoutput)
Results of Dickey Fuller Test: Test Statistic -6.625165e+00 p-value 5.901289e-09 #Lags Used 4.300000e+01 Number of Observations Used 1.659300e+04 Critical Value (1%) -3.430744e+00 Critical Value (5%) -2.861714e+00 Critical Value (10%) -2.566863e+00 dtype: float64
Onze data is stationair want p-value is 5,901289 x 10^-9 dit is veel kleiner dan een significantie level van 0.05 dus is de data stationair
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
# Seasonal ACF plot
plot_acf(df_time['cnt'], lags=40, ax=ax1)
ax1.set_title('Seasonal Autocorrelation Function (ACF)')
ax1.set_ylim(0.2, -0.2)
# Seasonal PACF plot
plot_pacf(df_time['cnt'], lags=40, ax=ax2)
ax2.set_title('Seasonal Partial Autocorrelation Function (PACF)')
ax2.set_ylim(0.2, -0.2)
plt.show()
#Het splitsen van de data. Dit gaat anders dan bij het splitsen van niet-tijdreeksgegevens.
split_point = len(df_train) - (3 * 7 * 24)
train, test = df_train[0:split_point], df_train[split_point:]
#Het model trainen
sarimax_model = SARIMAX(train['cnt'], order=(5, 1, 6), seasonal_order=(0, 1, 1, 24))
sarimax_model = sarimax_model.fit(disp=True, maxiter=1000)
c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
warnings.warn("Maximum Likelihood optimization failed to "
De order hebben we gekozen door veel te proberen en uiteindelijke de beste te kiezen. Hieronder laten we zien welke order en seasonal_order we hebben gebruikt.
We hebben voor een p van 5 gekozen want als je bij de ACF de bolletjes telt tot het eerste bolletje dat in de betrouwbaarheids ding (?) zit, zijn dat er 5. De q van 6 hebben we op dezelfde manier geteld maar dan bij de PACF plot. De d is 1 omdat de data stationair is.
Voor seasonal error wisten we dat de S 24 is want dat is onze seasonal period.
En voor de PDQ hebben we veel geprobeerd en uiteindelijk zijn wij op 0,1,1 gekomen.
sarimax_model.plot_diagnostics(figsize=(15,12))
plt.show()
Model diagnostics
De functie plot_diagnostics genereert vier grafieken om de prestaties van het model te beoordelen. Deze grafieken behandelen gestandaardiseerde residuen, het histogram met de schatting van de kansdichtheid (KDE), de normale kwantiel-kwantiel (Q-Q) plot en het correlogram.
Om het model als geschikt te beschouwen, moeten aan de volgende voorwaarden worden voldaan.
Gestandaardiseerde Residuen:
Histogram met KDE-schatting:
Normale Q-Q Plot:
Correlogram (ACF Plot):
(Artley, 2022)
In onze grafieken zie je het volgende:
standardized residual
Histogram plus estimated density
Normal Q-Q
Correlogram
forecast = sarimax_model.get_forecast(steps=len(test))
predictions = forecast.predicted_mean
confidence_intervals = forecast.conf_int()
lower_limits = confidence_intervals.iloc[:, 0]
upper_limits = confidence_intervals.iloc[:, 1]
# Plotting
plt.figure(figsize=(12, 6))
plt.plot(train.index, train['cnt'], label='Train')
plt.plot(test.index, test['cnt'], label='Test', color = 'g')
plt.plot(predictions.index, predictions, label='Predictions', color='r')
# Fill area between confidence intervals
plt.fill_between(predictions.index, lower_limits, upper_limits, color='pink', alpha=0.5)
plt.legend()
plt.title('SARIMA Predictions with 95% Confidence Intervals')
plt.ylim(-50, 900)
plt.xlim(pd.to_datetime('2012-11-07'), pd.to_datetime('2012-11-30'))
plt.show()
Aan de hand van de grafiek kunnen we zien dat het model niet perfect voorspelt maar wel in de goede richting zit. Want wanneer de test (groen) lijn omhoog gaag gaat de lijn van de predictions (rood) ook omhoog en hetzelfde geld voor wanneer de test lijn omlaag gaat de predictions lijn ook omlaag.
# Assuming you have a test dataset called 'test' with the corresponding 'cnt' values
test['predicted'] = sarimax_model.predict(start=test.index[0], end=test.index[-1], dynamic=True)
# Calculate the squared differences between predicted and actual values
squared_diff = (test['cnt'] - test['predicted'])**2
# Calculate the mean squared error
mse = np.mean(squared_diff)
# Calculate the root mean squared error
rmse = np.sqrt(mse)
print(f'Root Mean Squared Error (RMSE): {rmse}')
Root Mean Squared Error (RMSE): 128.71632798640323
C:\Users\crazy\AppData\Local\Temp\ipykernel_23884\2082032065.py:2: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy test['predicted'] = sarimax_model.predict(start=test.index[0], end=test.index[-1], dynamic=True)
Nu gaan we sarimax trainen op de gehele train dataset, zodat deze op de test dataset kan worden gebruikt.
# Model trainen op gehele train set
sarimax = SARIMAX(df_train['cnt'], order=(5, 1, 6), seasonal_order=(0, 1, 1, 24))
sarimax = sarimax.fit(disp=True, maxiter=1000)
Nu gebruiken we SARIMAX op de test data en maken we een csv bestand aan voor Kaggle.
# Haal de voorspellingen op
sar_pred = sarimax_model.get_forecast(steps=len(df_test))
sar_pred = sar_pred.predicted_mean
# Maak de voorspellingen tot lijst
sar_pred = sar_pred.tolist()
# Aanmaken df voor csv
test_predictions_df = pd.DataFrame(
{'date_hour': data_test['date_hour'],
'cnt': sar_pred})
# Aanmaken csv
test_predictions_df.to_csv(
f'Kaggle Submissions/vs_SARIMAX.csv',
index=False)
Prophet is een methode voor het voorspellen van tijdreeksgegevens. Het maakt gebruik van een statische model, wat betekent dat het de voorspelde uitkomst beschouwt als de som van verschillende componenten. Deze componenten omvatten niet-lineaire trends die worden aangepast aan jaarlijkse, wekelijkse en dagelijkse seizoensinvloeden, evenals speciale effecten rond vakantieperiodes. Prophet is vooral effectief bij tijdreeksen die sterke seizoenseffecten vertonen en waarbij historische gegevens beschikbaar zijn voor meerdere seizoenen. Een belangrijk kenmerk van Prophet is zijn robuustheid tegen ontbrekende gegevens en het vermogen om goed om te gaan met trendverschuivingen en uitschieters. (Prophet, z.d.)
De code die we hebben gebruikt komt uit de bron: (Quick start, 2023)
#kolom namen aanpasse naar datastamp (ds) en y (target value)
df1 = df.rename(columns={'date_hour': 'ds', 'cnt': 'y'})
# Opstellen van Prophet model
prophet_model = Prophet()
# model fitten
prophet_model.fit(df1)
# dataframe maken voor predictions
# periods is 24 want seasonal period is 24 uur
future = prophet_model.make_future_dataframe(periods=24)
# voorspelling maken
forecast = prophet_model.predict(future)
# voorspelling plotten
fig = prophet_model.plot(forecast)
plt.title('Prophet Forecast')
plt.show()
21:52:06 - cmdstanpy - INFO - Chain [1] start processing 21:52:08 - cmdstanpy - INFO - Chain [1] done processing c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:72: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result fcst_t = fcst['ds'].dt.to_pydatetime() c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:73: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result ax.plot(m.history['ds'].dt.to_pydatetime(), m.history['y'], 'k.',
Hier worden de componenten van de voorspelling laten zien op trend, week en dag niveau.
fig2 = prophet_model.plot_components(forecast)
plt.show()
c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:228: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result fcst_t = fcst['ds'].dt.to_pydatetime() c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:397: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result artists += ax.plot(df_y['ds'].dt.to_pydatetime(), seas[name], ls='-', c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:401: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result df_y['ds'].dt.to_pydatetime(), seas[name + '_lower'], c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:397: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result artists += ax.plot(df_y['ds'].dt.to_pydatetime(), seas[name], ls='-', c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:401: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result df_y['ds'].dt.to_pydatetime(), seas[name + '_lower'],
Dit is een interactive grafiek waar je per week, maand, half jaar of een jaar kan kijken naar de componenten van de voorspellingen
actual = df['cnt']
# 24 * 365 = 8760
predicted = forecast['yhat'].tail(8760).reset_index(drop=True)
# Bekijken van de RMSE
rmse = np.sqrt(mean_squared_error(actual[-8760:], predicted))
print(f'RMSE: {rmse}')
RMSE: 132.69287924706384
Nu maken we de voorspelling op de test set.
# Train dataset
df_train_reset = df_train.reset_index()
train_data = df_train_reset.rename(columns={'date_hour': 'ds', 'cnt': 'y'})
# Maken model
proph = Prophet()
proph.fit(train_data)
# Test dataset
test_data = df_test.reset_index()
prediction_data = test_data.rename(columns={'date_hour': 'ds'})
# Voorspellingen maken
forecast = proph.predict(prediction_data)
# Aanmaken df voor csv
test_predictions_df = pd.DataFrame(
{'date_hour': data_test['date_hour'],
'cnt': forecast['yhat']})
# Aanmaken csv
test_predictions_df.to_csv(
f'Kaggle Submissions/vs_Prophet.csv',
index=False)
22:13:19 - cmdstanpy - INFO - Chain [1] start processing 22:13:22 - cmdstanpy - INFO - Chain [1] done processing
Om alle modellen gemakkelijk te vergelijken op resultaat, is er een tabel gemaakt:
| Modellen | RMSE Score |
|---|---|
| Linear Regression: | 117.2269334691733 |
| Linear Regression Tuned: | 117.2269334691733 |
| KNN: | 53.38359939354527 |
| KNN Tuned: | 54.4182475724793 |
| Decision Tree: | 54.01557199126803 |
| Decision Tree Tuned: | 49.12543027929688 |
| Random Forest: | 37.84325840402451 |
| Random Forest Tuned: | 37.78426103437784 |
| XGBoost: | 37.458258453583944 |
| XGBoost Tuned: | 36.402051993205 |
| Ensemble, Stacking: | 40.11728088527348 |
| Ensemble, Stacking Tuned: | 37.7191700700007 |
| Hybrid Model: | 36.66069528728783 |
| SARIMA(X): | 128.71632798640323 |
| Prophet: | 132.69287924706384 |
Uit de bovenstaande tabel is te zien dat de laagste Root Mean Squared Error is voortgekomen uit de voorspellingen van het getunede XGBoost model. Het hybride model, die ook gebruik maakt van XGBoost, heeft de op een na laagste score op de test data. De parameters van het XGBoost waren als volgt:
xgb_tuned = XGBRegressor(n_estimators=300,
learning_rate=0.1,
max_depth=5,
random_state=SEED,
n_jobs=-1)
Op Kaggle is dit echter andersom en ook met een groter verschil tussen de scores. Dat deze modellen zo dicht bij elkaar liggen is niet raar, aangezien ze beide gebruik maken van het XGBoost model. Ook zijn de residuals een goede extra parameter om de voorspellingen te maken, aangezien ze data meegeven van een ander model. Hieronder zijn de scores voor op Kaggle te vinden.

Van de features, hebben de kolommen 'Uur' en 'is_weekend' door het hele proces de meeste invloed gehad op de score van de voorspellingen. Zeker in combinatie met elkaar, omdat de uren in de weekenden anders reageren dan tijdens de werkdagen. Daarnaast is het ook te zien dat de dagelijkse sinus beweging veel impact heeft. Dit komt door het gevonden dagelijkse seizoens patroon, die door middel van Fourier Analyse is verkregen.
In het EDA, met name in het time-series gedeelte, was dit ook duidelijk te zien. De grafiek over de verhuuraantallen ten opzichte van de uren en het weekend gaf een duidelijk beeld over de verschillen van deze periode en de potentie van deze kolommen. Wat ook opviel, is dat er in de top van de feature importances alleen zelf toegevoegde kolommen aanwezig waren.
Voor advies raden wij aan dat er meer wordt geadverteert op plekken waar het rondom de spits druk is. Daarnaast is het handig om ook op populaire locaties te adverteren. In het geheel zorgt dit ervoor dat meer mensen bekend worden met de service, wat leidt tot meer mogelijkheden tot verhuur. Daarnaast is het ook handig om dagelijks voor te bereiden op de drukte rond de juiste tijd. Omdat dit in de werkdagen verschilt van het weekend, is het handig om op de juiste tijden voorbereid te zijn.
We denken wel dat uiteindelijk een tijdseries model beter is om te voorspellen aangezien deze echt de toekomst voorspelt. Maar onze modellen doen het niet zo goed met de RMSE score. Vandaar dat wij hem nu nog niet aanbevelen maar als deze verbeterd is en hij doet het goed, raden wij zeker aan om een tijdseries model zoals bijvoorbeeld SARIMA(X) te gebruiken.
Voor alle onderwerpen is gebruik gemaakt van de Datacamps, de lesopgaven en de powerpoints.